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The topic of this presentation is the determination of orbits using Picard iteration, v/hich is 
a direct extension of the classical method of Picard that has been used in finding approxi- 
mate solutions of nonlinear differential equations for a variety of problems. The application 
of the Picard method of successive approximations to the initial value and the two-point 
boundary value problems is given below. 

The initial value problem, 

z=F(z, t) 

z(t,) = c, n = l,2, ... 

is solved '^y means of the iterative scheme, 

z„(t,) = c. 

whereby the nth approximation to the state vector z is computed from the (n - 1) approxi- 
mation (beginning with some initial approximation, z^ft)). The computation involves per- 
forming a simple integration or quadrature at each stage of the iterative process. The con- 
stant of integration is chosen such that each iterate satisfies the initial condition, z (t, ) = c. 

Similarly, the two-point boundary value problem, 

z = F (z, t) 

g(z (t,),z(tj)) = 0, 

is solved by means of the iterative scheme, 

z„ = F(Vi-0 

g(z„(t,),z^(tj)) = 0, n = l,2,... 
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where g denotes certain constraints or boundary conditions that must be satisfied by the 
solution at t, and tj . 

This paper presents an investigation of the suitability of this type of iterative scheme for tne 
problem of estimating orbits. In the estimation problem, we again have to solve the differ- 
ential equation for the state. However, we are now given a set of imprecise observations, 
y(tj ) j” , wfiich are, in general, nonlinear functions of the state. The problem is to find 
a z(t) which comes close (in some sense) to satisfying the observations. The estimation 
problem can thus be written as 

z = F (z, t) 

subject to the condition that z(t) minimizes some function Q, where 

Q = Q [z (t, ), z (tj) z (t„), y (t , ). }• (tj), . . . , y (t^)] . 

For instance, Q could be the sum of the squares of the residuals (the residual being the 
difference between the observed and ine computed values of the observations). Again, the 
iteration is set up, 

now choosing the constant of integration for z„(t) such that 

Q [z„ (t, ). (t j), (t„ ), y (t ,)...., y (t^ )1 

is minimized. As before, each stage of the iteration requires that a quadrature is performed. 
In order to perform the integration, we have chosen to approximate the right-hand side, 

F(Zj| , , t), by n series of Chebyshev polynomials, T^(t), and to integrate term-by-term. 

Thus, the approximation 

j=0 

is made where r = 1 - 2 - t, ) is normalized tme. The coefficients of the series are 

determined readily using the orthogonality properties from the relation 

k=0 

where tj = cos (jir/N). The integration can be performed by manipulating these coefficients, 
giving 

N , 
j=0 
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where 


aj=(bj., forj = l,2,...N 

and a,, = the arbitrary constant of integration that is determined from the requirement that 
Q be minimized. 

We will now consider the application of this technique to an orbit estimation prob:';;n., 
namely, the determination of the orbit of an earth satellite from observations o( the ange 
and range rate by several tracking stations. 

There are no external forces except the gravity field of the earth. The equations o' motion 
for the satellite are written as a set of six, first-order equations for the six state variables; 

X, y, z, u, V, and w. For instance, R,(t.) represents the range observed at t, from the fifth 
tracking station, and range rate is the time derivative of the range denoted by R,(t|): these 
are given. For the computer simulation presented here, we have taken the nominal values of 
the range and the range rate and added random noise. The standard deviations assumed are 
30 meters in range and about 50 cm per second in range rate. 

For this particular problem, the Q function should be considered. This is the sum of the 
weighted squares of the residuals; 

$ m 

-EE 

s=l 1=1 
*1 “l 

w., 

5=1 i=l 

where 

g,(tj) = |(x(t,)-x,(t.)l^ 

and 

x,(t,), yj(tj), Zj(t,) = the fifth station coordinatei a* t,, 

R,(tj) = the range from the particular tracking station, 

g, = the range computed using the piesent iterate, 

• 

Rj(tj) = observed range rate, and 

= weight for the particular observation. 
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We can see that is the function of the state at tj, and represents the coordinates of the 
station itself. The minimization of this function Q with respect to the integration constant 
is carried out using Newton's method. 

Table 1 provides the data for the three tracking stations we have assumed. The longitude 
and latitude, number of observations, and the total interval considered— about 1 500 sec- 
onds-are listed. The elliptic orbit considered has ar. ;centricity of 0.0557, with a semi- 
major axis of about 7178 km. This corresponds to a perigee at 400 km altitude, an apogee 
at 1200 km altitude, and an inclination of 20'* to the equatorial plane. 

Tab’s 1 

Station and Orbit Data 



Some of the results obtained are depicted in figure 1 , which shows the error as a function 
of time for different iterations. The initial guess is off by about 70 km from the true solu- 
tion, and with two iterations the error is brought down to something like 1 5 km. 

In order to see the convergence properly, the log of absolute error is plotted as a function 
of the iteration number in figure 2. The error is plotted for several points along the trajec- 
tory. The first point is T = 0, and the last point is T = 1 5 1 3. 1 . Jt is seen that the conver- 
gence is linear as was expected because of its relation to the classical Picard method. 

One of the disadvantages of this method is that the convergence is obtain;?d only for arcs 
of length less than ab it one-third of a revolution. Another problem is that of developing 
a procedure for con ng tv/o successive arcs. 
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Figure 1. Error as a function of time for different iterations 
(elliptic orbit). 


On the other hand, there are many advantages: For example, the method is quite simple and 
does not require the linear perturbation equations. The solution is in the form of a poly- 
nomial and is convenient to store. There is no interpolation required, should the solution 
be needed at some intermediate point. The error in representing this by a polynomial can, 
of course, be estimated by observing the last few terms of the polynomial. Also, the method 
is not sensitive to a poor initial guess. As we have seen here, the initial guess was off by 
about 70 km and even then the process will converge without difficulty. 

The question is, is this method reaUy more efficient than a more traditional approach to 
solving a problem? Actually, we have not compared it with the existing methods. However, 
we have some estimates on the time taken for this particular problem-it takes about 30 
seconds on the IBM 360/65. 


LOG OF ABSOLUTE ERROR 



Figure 2. Convergence of the method for an elliptic orbit. 


6 


A FINITE ELEMENT GRAVITY FIELD? 

J. L. Jiinkins 
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The paper details an approach for constructing a globally valid, piecewise continuous family 
of locall;' valid potential functions. The thesis put forth here is that the higher frequency 
terms of the geopotential can be more efficiently computed from such locally convergent 
approximations (typically three variable power series of order less than 5) than from any 
globally convergent gravity representation. This approach appears to be a step in the direc- 
tion of voiding the recent trend that the better model we have of the geopotential, the more 
expensive it is to integrate orbits with it! Numerical experiments conducted thus far confirm 
the validity of the approach and that acceleration errors of 0 (10''® km/s^ ) are achievable. 
The trade-ofis involved in selecting the finite element shape and size, and the order of the 
loca' approximations versus resulting accuracy, computational speed, and storage require- 
ment., are currently under study. 
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MEAN RATES OF THE ORBITAL ELEMENTS OF A SATELLITE 
PERTURBED BY A LENS SHAPED MASCON 


M. t. Ananda 
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A set of mean orbital rates are computed for a satellite perturbed by a lens shaped mascon. 
A disturbing potential in terms of the orbital elements of the satellite and the mascon 
parameters is developed. The partial derivatives of the potential with respect to the orbital 
elements are formed. These partials are averaged over the period of the satellite orbit to 
eliminate the short periodic terms. The averaged partials are substituted into the variation 
of parameters equations to give the mean orbital rates. In the limiting case, when the radius 
of the lens shaped mascon reaches zero the mean orbital rate due to a point mass is ob- 
tained. The orbital rates developed by the method described here are compared against the 
rates obtained by numerical differencing. The method developed here is used to reduce the 
Apollo-15 and -16 subsatellite data for lunar farside gravity determination. 
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